home *** CD-ROM | disk | FTP | other *** search
/ PC-Blue - MS DOS Public Domain Library / PC-Blue MS-DOS Public Domain Library - NYACC.iso / vol270 / contour.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-12-16  |  14.9 KB  |  538 lines
  1. /*    contour - find contours of a tabulated function of two variables 
  2.  
  3.     history...
  4.         28 Jun 86    Printing output file name.
  5.                     Optionally listing tiles visited.
  6.         26 Jun 86    User may specify a voltage and get estimated
  7.                     capacitance per unit length.
  8.  
  9.     bugs...
  10.         Needs better approximation for the true location of a contour
  11.         than a linear interpolation?
  12.  
  13.     method...
  14.         The rectangular grid is divided into triangles, or "tiles",
  15.         numbered as follows:
  16.  
  17.     (xmin,ymin)                                            (xmax,ymin)
  18.             -------------------------------        ----------
  19.             | 0  /| 2  /| 4  /| 6  /| 8  /|          /|2w-2/|
  20.             |   / |   / |   / |   / |   / |         / |   / |
  21.             |  /  |  /  |  /  |  /  |  /  |           |  /  |
  22.             | /   | /   | /   | /   | /   | /         | /   |
  23.             |/  1 |/  3 |/  5 |/  7 |/  9 |/          |/2w-1|
  24.             ---------------------------------         -------
  25.             | 2w /|2w+2/|    /|    /|    /|                /|
  26.             |   / |   / |   / |   / |   / |               / |
  27.             |  /  |  /  |  /  |  /  |  /                    |
  28.             | /   | /   | /   | /   | /  
  29.             |/2w+1|/    |/    |/  
  30.             -------------------
  31.             |    /|    /|              ...
  32.     
  33.                                     |/    |/    |/    |/    |
  34.                                ------------------------------
  35.                              /|    /|    /|    /|    /|    /|
  36.                             / |   / |   / |   / |   / |   / |
  37.                         |  /  |  /  |  /  |  /  |  /  |  /  |
  38.                           | /   | /   | /   | /   | /   | /   |
  39.                         |/    |/    |/    |/    |/    |2hw-1|
  40.                       -----------------------------------------
  41.                                                     (xmax,ymax)
  42.  
  43.         Each tile is assigned edge numbers as follows:
  44.  
  45.                       2        
  46.                    -------              /|          
  47.                   |      /             / |
  48.                   |     /             /  |
  49.                 1 |    /  0       0  /   | 1
  50.                   |   /             /    |
  51.                   |  /             /     |
  52.                   | /             /      |
  53.                   |/              -------
  54.                                      2
  55.                        
  56.         For each contour level, MAIN scans the grid to find a triangle
  57.         that the contour crosses, then calls FOLLOW.  FOLLOW follows
  58.         the contour, marking each tile it crosses, until it leaves the
  59.         grid or closes on itself.  For each tile crossed, FOLLOW calls
  60.         RECORD_POINT which calculates the point at which the contour
  61.         leaves the tile (using linear interpolation) and records it in
  62.         a stack.  MAIN prints the stack contents in reverse order and,
  63.         if the contour left the grid, calls FOLLOW again and prints the
  64.         points along the rest of the contour.
  65.  
  66.     timing...
  67.         this .bat file:
  68.             time <cr >>times
  69.             contour v.pot -n 3
  70.             time <cr >>times
  71.  
  72.         gave these timings:
  73.             Current time is 13:25:55.80
  74.  
  75.             Current time is 13:27:47.63 -> 1.8619
  76.             Current time is 13:36:45.47
  77.  
  78.             Current time is 13:38:25.56 -> 1.6691 changed tests in crosses()
  79.  
  80.             Current time is 14:11:03.38
  81.  
  82.             Current time is 14:12:43.15 -> 1.6602 eliminated a multiply in record_point()
  83.  
  84. */
  85.  
  86. #include <stdio.h>
  87. #include <math.h>
  88.  
  89. int next_edge[6]={2,0,1,1,2,0};
  90. int etest[6];
  91. int next_tile[3];
  92.  
  93. #define VERSION "0.2"
  94.  
  95. #define maximum(a,b) (a)>(b)?(a):(b)
  96. #define BUFSIZE 128
  97. #define LEFT_SIDE(t)   ((t)%w2==0)
  98. #define RIGHT_SIDE(t)  ((t)%w2==w2-3)
  99. #define TOP_SIDE(t)    (!((t)&1)&&(t<w2))
  100. #define BOTTOM_SIDE(t) (((t)&1) && ((t)>=bottom_row))
  101.  
  102. double adjust();
  103.  
  104. double xmin;
  105. double ymin;
  106. double xmax;
  107. double ymax;
  108. double xp,yp;            /* last points sent to output file */
  109. double zmin=1.e30;
  110. double zmax=-1.e30;
  111. double ex,ey;            /* grid spacing in x and y direction */
  112. double slope;            /* ey/ex */
  113. double *x;                /* pointer to data area for contour x values */
  114. double *y;                /* pointer to data area for contour y values */
  115. double *z;                /* pointer to data area for tabulated heights */
  116. double xscale,yscale,zscale;
  117. double aspect,correct,eps;
  118. double    z1;            /* minimum z contour value */
  119. double    z2;            /* maximum z contour value */
  120. double adj;            /* scale factor for printing contour value */
  121. double integral;    /* integral of gradient along a contour */
  122. double voltage;        /* voltage between boundaries given in -v switch */
  123.  
  124. int contours=6;        /* number of contour levels */
  125. int debugging=0;
  126. int calculating_capacitance=0; /* nonzero if user specified voltage */
  127. int continuing;        /* zero on 1st call to record_point() for each contour */
  128. int error;            /* code returned by fprintf() */
  129. int h;                /* height of grid in cells */
  130. int w;                /* width of grid in cells */
  131. int w2;                /* w * 2 */
  132. int hw;                /* height*width */
  133. int hw2;            /* height*width*2 */
  134. int bottom_row;        /* tile number for first tile in bottom row */
  135. int kind=0;            /* nonzero for log contours */
  136. int np=0;            /* number of contour points stored so far */
  137. int num_points=0;    /* number of contour points which can be stored */
  138. z_arguments=0;        /* # arguments of -z switch */
  139. char outname[40]="out.cnt";
  140. char *mark;            /* nonzero if this tile has been visited */
  141. char buf[BUFSIZE];    /* text I/O buffer */
  142. char format[50];    /* format for printing scale factor */
  143.  
  144. FILE file;
  145. FILE ofile;
  146.  
  147.  
  148. main(argc,argv) int argc; char **argv;
  149. {    int i,j,k,xx,yy,zz,xx2,yy2,rlz,t,e,t0,e0,more;
  150.     double yval,h0,h1,h2;
  151.     char *s;
  152.     double zval;
  153.  
  154.     read_data(argc,argv);
  155. /*
  156.     attributes();
  157.     show_z();
  158. */
  159.     etest[0]=0;
  160.     etest[1]=1;
  161.     etest[2]=w;
  162.     etest[3]=w+1;
  163.     etest[4]=w;
  164.     etest[5]=1;
  165.     next_tile[0]=1;
  166.     next_tile[1]=-1;
  167.     next_tile[2]=1-2*w;
  168.  
  169.     rlz=contours+2;
  170.     scale_one(zmin,zmax,&z1,&z2,rlz,&contours,kind);
  171. #ifdef xxx
  172.     printf("main: zmin=%f  zmax=%f  kind=%d \n", zmin, zmax, kind);
  173.     printf("main: z1=%f z2=%f contours=%d \n",z1,z2,contours);
  174. #endif
  175.  
  176.     if(!kind) adj=adjust(format,z1,z2,contours);
  177.     zval=z2;
  178.     for(i=0; i<=contours; i++,zval -= (z2-z1)/contours)
  179.         {if(!kind) sprintf(buf,format,zval*adj);
  180.         else sprintf(buf," 1e%1.0f",zval);
  181.         printf("contour level %s",buf);
  182. /*
  183.         if(zval<=zmin || zval>=zmax) continue;
  184. */
  185.         for (t=0; t<hw2; t++) mark[t]=0; /* clear all marks */
  186.         for (t=0; t<hw2-w2; t++)
  187.             {integral=0.;
  188.             if(t%w2==w2-2)
  189.                 {t++;
  190.                 continue;  /* skip phantom tiles along right edge */
  191.                 }
  192.             if(!mark[t] && crosses(t,zval))
  193.                 {if(t&1) h0=z[t/2+w+1]; else h0=z[t/2];
  194.                 h1=z[t/2+1]; h2=z[t/2+w];
  195.                 if(h0<zval && zval<=h1) e=2;
  196.                 else if (h1<zval && zval<=h2) e=0;
  197.                 else e=1;
  198.                 t0=t;
  199.                 continuing=0;
  200. /* printf("main: picking up %g contour at tile %d, edge %d\n",zval,t,e); */
  201.                 more=follow(t,e,1,zval);
  202.                 if(more)
  203.                     {for (j=np; j; )
  204.                         {j--;
  205.                         error=fprintf(ofile,"%8.3g %8.3g\n",x[j],y[j]);
  206.                         if(error==-1) quit();
  207.                         }
  208.                     xp=x[0]; yp=y[0];
  209.                  /* we'll now be going opposite direction along contour... */
  210.                     integral= -integral;
  211.                     np=0;
  212.                     if(h1<zval && zval<=h0) e=2;
  213.                     else if (h2<zval && zval<=h1) e=0;
  214.                     else e=1;
  215. /* printf("main: restarting %g contour at tile %d, edge %d\n",zval,t0,e); */
  216.                     follow(t0,e,0,zval);
  217.                     integral= -integral;
  218.                     flush_contour(zval);
  219.                     }
  220.                 else flush_contour(zval);
  221.                 }
  222.             }    /* all cells checked for this contour level */
  223.         }
  224.     fclose(ofile);
  225. }
  226.  
  227. follow(t,e,high_on_left,zval) int t,e,high_on_left; double zval;
  228. {    int test;
  229.     while(1)
  230.         {if(debugging)
  231.             printf("contour entering tile %3d, edge %d %11s with high ground on %s\n",
  232.             t,e,
  233.             (e==0)?"(DIAGONAL),"
  234.                 :((t&1)?
  235.                     ((e==1)?"(LEFT),":"(BOTTOM),")
  236.                     :((e==1)?"(RIGHT),":"(TOP),")
  237.                 ),                        
  238.             high_on_left?"LEFT":"RIGHT");
  239.         if(t&1)                            /*    lower right triangle */
  240.             test=etest[e+3];
  241.         else                            /*    upper left triangle */
  242.             test=etest[e];
  243.         if(high_on_left ^ (z[t/2+test]>zval)) e=next_edge[e+3];
  244.         else                                  e=next_edge[e];
  245.         record_point(t,e,zval);
  246.         if((e==1 && (LEFT_SIDE(t) || RIGHT_SIDE(t))) ||
  247.             (e==2 && (TOP_SIDE(t) || BOTTOM_SIDE(t)))   ) return 1;
  248.         if(t&1) t -= next_tile[e];
  249.         else    t += next_tile[e];
  250.         if(mark[t]) return 0;
  251.         mark[t]++;
  252.         }
  253. }
  254.  
  255. attributes()
  256. {    printf("tile attributes...\n");
  257.     int j,k,m;
  258.     char type[9];
  259.     m=0;
  260.     for (k=0; k<h-1; k++)
  261.         {for (j=0; j<w2-2; j++)
  262.             {type[0]=0;
  263.             if(LEFT_SIDE(m)) strcat(type,"L");
  264.             if (TOP_SIDE(m)) strcat(type,"T");
  265.             if (RIGHT_SIDE(m)) strcat(type,"R");
  266.             if (BOTTOM_SIDE(m)) strcat(type,"B");
  267.             printf("%2s ",type);
  268.             m++;
  269.             }
  270.         printf("\n");
  271.         m += 2;
  272.         }
  273. }
  274.  
  275. show_z()
  276. {    int i,j,k=0;
  277.     printf("h=%d  w=%d\n",h,w);
  278.     for (i=0; i<h; i++)
  279.         {for (j=0; j<w; j++)
  280.             printf("%5.2f ",z[k++]);
  281.         printf("\n");
  282.         }
  283. }
  284.  
  285. record_point(t,e,zval) int t,e; double zval;
  286. {    double xt,yt,xc,yc,d,f=0.,dist,dx,dy;
  287.     int t5,t5w;
  288.     t5=t/2;
  289.     t5w=t5+w;
  290. /*    printf("record_point: tile %d  edge %d  zval %8.2g",t,e,zval); */
  291.     if(np==num_points)
  292.         {xt=x[np-1]; yt=y[np-1];
  293.         flush_contour(z);
  294.         x[0]=xt; y[0]=yt; np++;
  295.         }
  296.  
  297.     if(t&1) /* lower triangle */
  298.         {switch(e)
  299.             {case 2:
  300.                 d=z[t5w]-z[t5w+1];
  301.                 dx=-d*slope; 
  302.                 dy=(z[t5w+1]-z[t5+1]);
  303.                 if(d) f=(z[t5w]-zval)/d;
  304.                 xc=t5%w+f; yc=t5/w+1;
  305.                 break;
  306.             case 1:
  307.                 d=z[t5+1]-z[t5w+1];
  308.                 dx=(z[t5w+1]-z[t5w])*slope;
  309.                 dy=-d;
  310.                 if(d) f=(z[t5+1]-zval)/d;
  311.                 xc=t5%w+1; yc=t5/w+f;
  312.                 break;
  313.             case 0:
  314.                 d=z[t5+1]-z[t5w];
  315.                 dx=(z[t5w+1]-z[t5w])*slope;
  316.                 dy=(z[t5w+1]-z[t5+1]);
  317.                 if(d) f=(z[t5+1]-zval)/d;
  318.                 xc=t5%w+1.-f; yc=t5/w+f;
  319.                 break;
  320.             }
  321.         }
  322.     else /* upper triangle */
  323.         {switch(e)
  324.             {case 2:
  325.                 d=z[t5]-z[t5+1];
  326.                 dx=-d*slope;
  327.                 dy=(z[t5w]-z[t5]);
  328.                 if(d) f=(z[t5]-zval)/d;
  329.                 xc=t5%w+f; yc=t5/w;
  330.                 break;
  331.             case 1:
  332.                 d=z[t5]-z[t5w];
  333.                 dx=(z[t5+1]-z[t5])*slope;
  334.                 dy=-d;
  335.                 if(d) f=(z[t5]-zval)/d;
  336.                 xc=t5%w; yc=t5/w+f;
  337.                 break;
  338.             case 0:
  339.                 d=z[t5+1]-z[t5w];
  340.                 dx=(z[t5+1]-z[t5])*slope;
  341.                 dy=(z[t5w]-z[t5]);
  342.                 if(d) f=(z[t5+1]-zval)/d;
  343.                 xc=t5%w+1.-f; yc=t5/w+f;
  344.                 break;
  345.             }
  346.         }
  347.     x[np]=xt=xc*ex+xmin;
  348.     y[np]=yt=yc*ey+ymin;
  349.     if(continuing)
  350.         integral += dx*(yt-yp)-dy*(xt-xp);
  351.     else
  352.         continuing=1;
  353. /*    dx=xp-xt; dy=yp-yt; dist=dx*dx+dy*dy;
  354.     if(dist>eps) */
  355.         np++;
  356.     xp=xt; yp=yt;
  357. /*    printf(" x,y=(%7.2g,%7.2g)\n",x[np],y[np]); */
  358. }
  359.  
  360. crosses(t,zval) int t; double zval;
  361. {    int gt=0, le=0;
  362.     if(z[t/2+1]<zval) gt=1; else le=1;
  363.     if(z[t/2+w]<zval) gt=1; else le=1;
  364.     if(t&1)
  365.         {if(z[t/2+w+1]<zval) gt=1; else le=1;
  366.         }
  367.     else
  368.         {if(z2=z[t/2]<zval) gt=1; else le=1;
  369.         }
  370.     return (gt&le);
  371. /*
  372.     int does_cross;
  373.     does_cross=(gt&le);
  374.     printf("crosses: contour at %g %s cross tile %d (%4.2f, %4.2f, %4.2f)\n",
  375.         zval,does_cross?"DOES":"DOES NOT",t,z[t/2+1],z[t/2+w],z2);
  376.     return does_cross;
  377. */
  378. }
  379.  
  380. read_data(argc,argv) int argc; char **argv;
  381. {    int i,j,k,length;
  382.     double xx,yy,d,*pd,sum;
  383.     char *s,*t;
  384.  
  385.     if(argc>1 && streq(argv[1],"?")) help();
  386.     if(argc<=1 || *argv[1]=='-') file=stdin;
  387.     else
  388.         {if(argc>1)
  389.             {argc--; argv++;    /* advance to file name */
  390.             file=fopen(*argv,"r");
  391.             if(file==0) {printf("file %s not found\n",*argv); exit();}
  392.             strcpy(outname,*argv);
  393.             if(s=index(outname,'.')) *s=0;
  394.             strcat(outname,".cnt");
  395.             }
  396.         else help();
  397.         }
  398.     argc--; argv++;        /* advance to first switch if any */
  399.     while(argc>0)
  400.         {i=get_parameter(argc,argv);
  401.         argv=argv+i; argc=argc-i;
  402.         }
  403.     ofile=fopen(outname,"w"); 
  404.     if(ofile==0) {printf("can\'t create output file\n"); exit(1);}
  405. /*
  406.     sample data file header...
  407. 0.  3.  10      minimum and maximum y values, width of grid
  408. 0.  3.  10      minimum and maximum y values, height of grid
  409.  
  410. */
  411.     if(mygets(buf,BUFSIZE,file,ofile)==0)
  412.             {printf("failure on reading input file header");
  413.             exit(1);
  414.             }
  415. /*    printf("1st input line:%s",buf); */
  416.     sscanf(buf,"%F %F %d",&xmin,&xmax,&w);
  417.     if(mygets(buf,BUFSIZE,file,ofile)==0)
  418.             {printf("failure on reading input file header");
  419.             exit(1);
  420.             }
  421. /*    printf("2nd input line:%s",buf); */
  422.     sscanf(buf,"%F %F %d",&ymin,&ymax,&h);
  423.     if(w<2 || h<2)
  424.         {printf("invalid data header... width or height too small\n");
  425.         exit(1);
  426.         }
  427.     ex=(xmax-xmin)/(w-1);
  428.     ey=(ymax-ymin)/(h-1);
  429.     if(ex==0. || ey==0.)
  430.         {printf("invalid data header... xmin==xmax or ymin==ymax\n");
  431.         exit(1);
  432.         }
  433.     slope=ey/ex;
  434.     if(ey<ex) eps=ey;
  435.     else eps=ex;
  436.     eps=eps*eps*1.e-4;
  437. /*
  438.     printf("xmin=%g  xmax=%g  ex=%g  w=%d\n",
  439.         xmin,xmax,ex,w);
  440.     printf("ymin=%g  ymax=%g  ey=%g  h=%d\n",
  441.         ymin,ymax,ey,h);
  442. */
  443.     if(h<2||w<2)
  444.         {fprintf(stderr,"parameter too small: height=%d, width=%d",h,w);
  445.         exit(1);
  446.         }
  447.     hw=h*w;
  448.     w2=w*2;
  449.     hw2=hw*2;
  450.     bottom_row=hw2-w*4;
  451.     num_points=4*(w+h);
  452.     x=malloc(num_points*sizeof(double));
  453.     y=malloc(num_points*sizeof(double));
  454.     z=malloc(hw*sizeof(double));
  455.     mark=malloc(hw*2);
  456.     if(x==0 || y==0 || z==0 || mark==0)
  457.         {fprintf(stderr,"can\'t allocate buffers"); exit();
  458.         }
  459.     k=0;
  460.     while(k<hw)
  461.         {if(fscanf(file,"%F",&z[k])==-1) exit();
  462.         switch(z_arguments)
  463.             {case 0:if(z[k]<zmin) zmin=z[k];
  464.             case 1:    if(z[k]>zmax) zmax=z[k];
  465.             case 2:    k++;
  466.             }
  467.         }
  468.     printf("writing output to %s\n",outname);
  469. }
  470.  
  471. /* get_parameter - process one command line option
  472.         (return # parameters used) */
  473. get_parameter(argc,argv) int argc; char **argv;
  474. {    int i;
  475.     if(streq(*argv,"-d")) {debugging=1; return 1;}
  476.     else if(streq(*argv,"-n"))
  477.         {if((argc>1) && numeric(argv[1])) contours=atoi(argv[1]);
  478.         return 2;
  479.         }
  480.     else if(streq(*argv,"-o"))
  481.         {argv++;
  482.         if(argc>1 && **argv!='-')
  483.             {strcpy(outname,*argv);
  484.             return 2;
  485.             }
  486.         fprintf(stderr,"missing argument for -o switch");
  487.         }
  488.     else if(streq(*argv,"-v"))
  489.         {i=get_double(argc,argv,1,&voltage,&voltage,&voltage);
  490.         if(i && fabs(voltage)>1.e-6) calculating_capacitance++;
  491.         return i;
  492.         }
  493.     else if(streq(*argv,"-z"))
  494.         {i=get_double(argc,argv,2,&zmin,&zmax,&zmax);
  495.         z_arguments=i-1;
  496.         if(z_arguments>2) z_arguments=2;
  497.         else if(z_arguments<0) z_arguments=0;
  498.         return i;
  499.         }
  500.     else gripe(argv);
  501. }
  502.  
  503. help()
  504. {    fprintf(stderr,"contour   version %s",VERSION);
  505.     fprintf(stderr,"\nusage: contour  [file]  [options]\n");
  506.     fprintf(stderr,"options are:\n");
  507.     fprintf(stderr,"  -n  num            approximate number of contours (default %d) \n",contours);
  508.     fprintf(stderr,"  -o  name           specify name of output file \n");
  509.     fprintf(stderr,"  -v  voltage        calculate capacitance assuming this \n");
  510.     fprintf(stderr,"                     voltage between the two boundaries\n");
  511.     fprintf(stderr,"  -z  [min [max]]    specify range of contours\n");
  512.     exit();
  513. }
  514.  
  515. flush_contour(zval) double zval;
  516. {    int i;
  517.     i=0;
  518.     while(1)
  519.         {error=fprintf(ofile,"%8.3g %8.3g",x[i],y[i]);
  520.         if(error==-1) quit();
  521.         if(++i>=np) break;
  522.         error=fprintf(ofile,"\n");
  523.         if(error==-1) quit();
  524.         }
  525.     printf("   Integral of (gradient dot normal) is %7.2g\n",integral/ey);
  526.     strcat(buf,"\n");
  527.     fputs(buf,ofile);
  528.     error=fprintf(ofile,"; integral of (gradient dot normal) was %7.2g",integral/ey);
  529.     if(error==-1) quit();
  530.     if(calculating_capacitance)
  531.         printf(" -> C =%7.2g pF/m\n",1e12*8.85e-12/voltage*integral/ey);
  532.     else
  533.         printf("\n");
  534.     integral=0.;
  535.     np=0;
  536. }
  537.  
  538. quit()
  539. {    fprintf(stderr,"output file error\n");
  540.     fclose(ofile);
  541.     exit(1);
  542. }
  543.